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Abstract. We determine the probabihty that an ultrahigh energy (above 5 • 10^® eV) 
proton created at a distance r with energy E arrives at earth above a threshold Ec- 
The clustering of ultrahigh energy cosmic rays suggests that they might be emitted by 
compact sources. We present a statistical analysis on the source density based on the 
multiplicities. The ultrahigh energy cosmic ray spectrum is consistent with the decay 
of GUT scale particles. By using a maximum likelihood analysis we determine the mass 
of these GUT scale particles. We consider the possibility that a large fraction of the 
ultrahigh energy cosmic rays are decay products of Z bosons which were produced in 
the scattering of ultrahigh energy cosmic neutrinos on cosmological relic neutrinos. 
Based on this scenario we determine the required mass of the heaviest relic neutrino as 
well as the necessary ultrahigh energy cosmic neutrino flux via a maximum likelihood 
analysis. 

1 Introduction 

The interaction of protons with the microwave background predicts a drop in 
the cosmic ray flux above the GZK cutoff wS-lO'^^ eV. The data shows no 
such drop. About 20 events even above 10^" eV were observed by a number of 
experiments. Since above the GZK energy the attenuation length of particles is a 
few tens of megaparsecs if an ultrahigh energy cosmic ray (UHECR) is observed 
on earth it was most probably produced in our vicinity. 

Section 2 studies the propagation and determines the probability /'(r, E, Ec) 
that protons created at distance r with energy E reach earth above a threshold 
Ec- Using this P one can give the observed spectrum by one numerical integration 
for any injection spectrum. 

It is an interesting phenomenon that the UHECR events are clustered. Usu- 
ally it is assumed that at these high energies the galactic and extragalactic 
magnetic fields do not affect the orbit of the cosmic rays, thus they should point 
back to their origin within a few degrees. In contrast to the low energy cosmic 
rays one can use UHECRs for point-source search astronomy. Recently, a sta- 
tistical analysis based on the multiplicities of the clustered events estimated 
the source density. In Section 3 we extend the above analysis. Our analytical 
approach gives the event clustering probabilities for any space, intensity and en- 
ergy distribution of the sources by using a single additional propagation function 
P{r,E;E,). 
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In Section 4 we study the scenario that the UHECRs are coming from decay- 
ing superheavy particles (SP) and we determine their masses mx by an analysis 
of the observed UHECR spectrum. 

Ultrahigh energy neutrinos (UHEz/) scatter on relic neutrinos (Riy) producing 
Z bosons, which can decay hadronically (Z-burst) In Section 5 we compare 
the predicted proton spectrum with observations and determine the mass of the 
heaviest Hv via a maximum likelihood analysis. 

The details of the presented results and a more complete reference list can 
be found in [|j|J|. 

2 Propagation of UHECR protons 

Using pion production as the dominant effect of energy loss for protons at en- 
ergies >10^^ eV, ref. calculated P{r, E, Ec) for three threshold energies. We 
extended the results of . The inelasticity of Bethe-Heitler pair production is 
small (w 10"'^), thus we used a continuous energy loss approximation for this 
process. The inelasticity of pion-photoproduction is larger (w 0.2 — 0.5) in the 
energy range of interest, thus there are only a few tens of such interactions during 
the propagation. Due to the Poisson statistics and the spread of the inelasticity, 
we will see a spread in the energy spectrum even if the injected spectrum is 
mono-energetic . 

In our simulation protons are propagated in small steps (10 kpc), and after 
each step the energy losses due to pair production, pion production and the 
adiabatic expansion are calculated. During the simulation we keep track of the 
current energy of the proton and its total displacement. We used the follow- 
ing type of parametrization P{r, E, Ec) = exp [—a • (r/1 Mpc)^] . Fig. |^ shows 
a{E/Ec) and b{E/Ec) for a range of three orders of magnitude and for five differ- 
ent Ec- Just using the functions oi a{E / E^) and b{E/Ec), thus a parametrization 
of P{r^ E, Ec), one can obtain the observed energy spectrum for any injection 
spectrum without additional Monte-Carlo simulation. 

Since P{r,Ep;E) is of universal usage, we have decided to make the latest 
numerical data for the probability distribution {—)dP{r, Ep-, E)/dE available for 
the public via the World- Wide- Web URL 

http:/ /www. desy.de/~uheci . 

The propagation function can be similarly determined for photons, though 
the necessary CPU power is approximately 300 times more than for protons. 
Therefore, we used the stochastic method to test a few cases. Usually the con- 
tinuous energy loss approximation was used. In this approximation, the energy 
(and number) of the detected photons is a unique function of the initial energy 
and distance, and statistical fluctuations are neglected. The processes that are 
taken into account are pair production on the diffuse extragalactic photon back- 
ground, double pair production and inverse Compton scattering of the produced 
pairs. The energy attenuation length of the photons due to these processes is 
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strongly influenced by the poorely known universal radio and infrared back- 
grounds. These uncertainties influences some of our results. 

A full simulation of the photon propagation function with all the statistical 
fluctuations will be the subject of a later work. 



3 Density of sources 

The arrival directions of the UHECRs measured by experiments show some 
peculiar clustering: some events are grouped within ^ 3°, the typical angular 
resolution of an experiment. Above 4 • 10^^ eV 92 cosmic ray events were de- 
tected, including 7 doublets and 2 triplets. Above 10^° eV, one doublet out of 14 
events were found ||]. The chance probability of such a clustering from uniform 
distribution is rather small |^,||. 

The clustered features of the events initiated an interesting statistical analysis 
assuming compact UHECR sources The authors found a large number, ~ 400 
for the number of sources within the GZK sphere. We generalize their analysis. 
The most probable value for the source density is really large; however, the 
statistical significance of this result is rather weak. 

The number of UHECRs emitted by a source of luminosity A during a period 
T follows the Poisson distribution. However, not all emitted UHECRs will be 
detected. They might loose their energy during propagation or can simply go 
to the wrong direction. For UHECRs the energy loss is dominated by the pion 
production in interaction with the cosmic microwave background radiation. In 
the previous section the probability function P(r, Ec) was calculated. 

The features of the Poisson distribution enforce us to take into account the 
fact that the sky is not isotropically observed. 

The probability of detecting k events from a source at distance r with energy 
E can be obtained by simply including the factor P{r, E, Ec)Ari/{ATTr'^) in the 
Poisson distribution: 

e^p[-~P{r,E,Ec)rij/r^] 

[P{r,E,E,)r^j/r^]\ (1) 

where we introduced j = XT A/ {Air) and Arj/ (47rr^), which is the probability that 
an emitted UHECR points to a detector of area A. The factor rj represents the 
visibility of the source, which was determined by spherical astronomy. We denote 
the space, energy and luminosity distributions of the sources by /o(x), c{E) and 
h{j), respectively. The probability of detecting k events above the threshold E^, 
from a single source randomly positioned within a sphere of radius R is 

Pfc = / dV p(x) / dE c{E) / dj h{j) X 
JSr J Jo 

e^l>[-Pir,E,Ec)Tjj/r^] .^x 
[P{r,E,Ec)vj/r \ . (2) 



Clustering, GUT scale and neutrino masses in Ultrahigh energy cosmic rays 



5 



Denote the total number of sources within the sphere of sufficiently large 
radius (e.g. several times the GZK radius) by N and the number of sources 
that gave k detected events by Nk- Clearly, N = Ni and the total number 
of detected events is A^e = X^o"*-^*- ^^^^ probability that for N sources the 
number of different detected multiplets are Nk is: 

P{N,{Nk})^ml[—P^^. (3) 

fc=o 

For a given set of unclustered and clustered events (iVi and N2,N3,...) invert- 
ing the P{N, {Nk}) distribution function gives the most probable value for the 
number 

Note, that Pk and then P(N, {Nk}) are easily determined by a well behaved 
four-dimensional numerical integration for any c{E), h(j) and p{r) distribution 
functions. In order to illustrate the uncertainties and sensitivities of the results 
we used a few different choices for these distribution functions. 

For c{E) we studied three possibilities. The most straightforward choice is 
the extrapolation of the 'conventional high energy component' cx E~^. Another 
possibility is to use a stronger fall-off of the spectrum at energies just below the 
GZK cutoff, e.g. oc E^^. The third possibility is to assume that UHECRs are 
some decay products of metastable superheavy particles [|lO 11 l^ , |l3| , p^ 15, l^ , |l7| 



or topological defects |18(] . The superheavy particles decay into quarks and gluons 
which initiate multi-hadron cascades through gluon bremstrahlung |]l^ , po|j2l] , |2^ , ^ . 

In the recent analysis |^ the authors have shown that for a fixed set of 
multiplets the minimal density of sources can be obtained by assuming a delta- 
function distribution for h{j). We studied both this limiting luminosity, h{j) = 
6{j — J*), and a more realistic one with Schechter's luminosity function, which 
can be given as: h{j)dj = h ■ {j / j*)'^''^^ exp{-j / j^)d{j / j^). 

The space distribution of sources can be given based on some particular 
survey of the distribution of nearby galaxies or on a correlation length vq char- 
acterizing the clustering features of sources. For simplicity the present analysis 
deals with a homogeneous distribution of sources. 

In order to determine the confidence level (CL) regions for the source densities 
we used the frequentist method [^. We wish to set limits on S, the source 
density. Using our Monte-Carlo based P{r, E, Ec) functions and our analytical 
technique we determined p{Ni, N2, N^, which gives the probability of 

observing Ni singlet, N2 doublet, N3 triplet etc. events if the true value of 
the density is S and the central value of luminosity is j». For a given set of 
{Ni,i = 1,2,...} the above probability distribution as a function of S and 
determines the 68% and 95% confidence level regions in the S — plane. 

Fig. 1^ shows the confidence level regions for one of our models (with in- 
jected energy distribution c{E) cx E~^; and Schechter's luminosity distribution: 
h(,3)dj oc (j/i*)~^'^^ exp(— jyj*)d(j'/j*)). The regions are deformed, thin ellipse- 
like objects. For this model our final answer for the density is 180^^gg|'jy^|^'' • 
10^'^ Mpc^^, where the first errors indicate the 68%, the second ones in the 
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Fig. 2. The la (68%) and 2a (95%) confidence level regions for j« and the source 
density (14 UHECR with one doublet). 

parenthesis the 95% CLs, respectively. The choice of Q -h(j)cx SQ)- and, e.g. 
E~'^ energy distribution gives much smaller value: 2.77^2^53(2^70)10"'^ Mpc~'^, 
which is in a quite good agreement with the result of Ref. ^ . 

4 Decay of GUT scale particles 

An interesting idea discussed by refs. ]l5| , p^|jr^ is that SPs could be the source 
of UHECRs. Note, that any analysis of SP decay covers a much broader class 
of possible sources. Several non-conventional UHECR sources produce the same 
UHECR spectra as decaying SPs. We study the scenario that the UHECRs are 
coming from decaying SPs and we determine the mass of this X particle mx by 
a detailed analysis of the observed UHECR spectrum. 

The hadronic decay of SPs yields protons. They are characterized by the frag- 
mentation function (FF) D{x, Q^) which gives the number of produced protons 
with momentum fraction x at energy scale Q. For the proton's FF at present 
accelerator energies we use ref. . We evolve the FFs in ordinary and in super- 
symmetric QCD to the energies of the SPs. This result can be combined with the 
prediction of the MLLA technique, which gives the initial spectrum of UHECRs 
at the energy mx (cf . Fig. . ^) . Similar results are obtained by . 

Depending on the location of the source -halo or extragalactic (EG)- and the 
model -SM or MSSM- we study four different scenarios. In the EG case protons 
loose some fraction of their energies, described by P{r, E, Ec). We compare the 
predicted and the observed spectrums by a maximum likelihood analysis. This 
analysis gives the mass of the SP and the error on it. 

Fig. U shows the measured UHECR spectrum and the best fit, which is ob- 
tained in the EG-MSSM scenario. 
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X 



Fig. 3. The quark FFs at Q^IO^*^ GoV for proton/pion in SM (solid/dotted lino) and 
in MSSM (dashed/dashed-dotted line). We change from logarithmic scale to linear at 
X = 0.01. 




Log (E [eV]) 

Fig. 4. UHECR data with their error bars and the best fit from a decaying SP. There 
are no events above 3 x 10^" eV (shown by an arrow). Zero event does not mean zero 
flux, but an upper bound for the flux. Thus, the experimental flux is in the "hatched" 
region with 68% GL. 
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To determine the most probable value for the mass of the SP we studied 4 
scenarios. Fig. || contains the Xmm values and the most probable masses with 
their errors for these scenarios. 

The UHECR data favors the EG-MSSM scenario. The goodnesses of the fits 
for the halo models are far worse. 
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Fig. 5. The most probable values for the mass of the decaying ultra heavy dark matter 
with their error bars and the total values. Note that 21 bins contain nonzero number 
of events and the fit has 3 free parameters. 



The SM and MSSM cases do not differ significantly. The most important 
message is that the masses of the best fits (EG cases) are compatible within the 
error bars with the MSSM gauge coupling unification GUT scale: mx = lO*" GeV, 



where b = U.Gti^. 



5 Z-burst scenario 



Already in the early eighties there were some discussions about the possibil- 
ity that the ultrahigh energy neutrino spectrum could have absorption dips at 
energies around El°^=M'^/ {2mi,^)=A.2 ■ 10^^ (1 eV/m^.) eV due to resonant 
annihilation with relic neutrinos of mass mi^, predicted by the hot Big Bang, 
into Z bosons of mass Mz [|2^,^. Recently it was realized that the same an- 
nihilation mechanism gives a possible solution to the GZK problem ||^. It was 
argued that the UHECRs above the GZK cutoff are from these Z-bursts. The 
Z-burst hypothesis for the ultrahigh energy cosmic rays was discussed in many 
papers |o|l|,||,|o|||,|n|,||,||3|§,0|8l||. 
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We compare this scenario with observations. 

The density distribution of Rt/s as hot DM follows the total mass distribu- 
tion; however, it is less clustered. Thus we, as opposed to practically all previous 
authors ||^,^,^, do not follow the unnatural assumption of having a relative 
overdensity of 10^ 10^ in our neighborhood (for an approach with lepton asym- 



metry see 32 ). Our quantitative results turned out to be rather insensitive to 
the variations of the overdensities within the considered range. 

We give the energy distribution of the produced particles in our lab system, 
which is obtained by Lorentz transforming the CM collider results. We included 
in our analysis the protons, which are directly produced in the Z-burst and 
appear as decay products of the neutrons. Photons were also taken into account. 
They are produced in hadronic Z decays via fragmentation into neutral pions, 
Z^7t'^ + X^2j + X. Electrons (and positrons) from hadronic Z decay are 
also relevant for the development of electromagnetic cascades. They stem from 
decays of secondary charged pions, Z^ir'^+X^e^+X. 

The UHECR flux from Z-bursts is proportional to the differential fluxes F^. 
of ultrahigh energy cosmic neutrinos. Unfortunately, the value of these fluxes 
is essentially unknown. In this situation of insufficient knowledge, we take the 
following approach concerning the flux of ultrahigh energy cosmic neutrinos of 
type i, F^. {E^,. , r). It is assumed to have the form 

F,^ {E,^ , r) = F,, {E,^ , 0) (1 + z)" , (4) 

where z is the redshift and where a characterizes the cosmological source evolu- 
tion. Note, however, that, independently of the production mechanism, neutrino 
oscillations result in a uniform mixture for the different mass eigenstates. 

The next ingredient of our analysis is the propagation of the protons and 
photons from cosmological distances. This propagation was described by the 
appropriate P{r, Ep, E) probability functions (see Section 2). 

Finally, we compare the predicted and observed spectrum and extract the 
mass of the relic and the necessary UHEj/ flux by a maximum likelihood 
analysis. Qualitatively, our analysis can be understood as follows. In the Z-burst 
scenario small relic ly mass needs large E^'^^ to produce a Z. Large El°^ results in 
a large Lorentz boost, thus large proton energy. In this way the detected energy 
determines the mass of the relic ly. The analysis is completely analogous to that 
of the previous section. The observed flux is a sum of two terms, namely the flux 
from Z-bursts and a conventional part with power-law behaviour in the energy. 
This power-law part might be produced in our galaxy (halo model) or it might 
be produced extragalactically (EG model). 

The Z-burst determination of the neutrino mass seems reasonably robust. 
Fig. ^ shows the summary of our relic neutrino mass determination. For a wide 
range of cosmological source evolution (a = —3 3), Hubble parameters h = 
0.61 0.9, f2M, ^A, -Zmax = 2 -f- 5, for variations of the possible relic neutrino 
overdensity in our GZK zone and for different assumptions about the diffuse 
extragalactic photon background, the results remain within the above error bars. 
The main uncertainties concerning the central values originate from the different 
assumptions about the background of ordinary cosmic rays. 
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In the case that the ordmary cosmic rays above 10^^'^ eV are protons and 
originate from a region within the GZK zone of about 50 Mpc ("halo"), the 
required mass of the heaviest neutrino seems to he between 2.1 eV< rrii, <6.7 eV 
at the 68% C.L. (a < 0), if we take into account the variations between the 
minimal and moderate universal radio background cases and the strong UHE7 
attenuation case. 

The much more plausible assumption that the ordinary cosmic rays above 
10^^'^ eV are protons of extragalactic origin leads to a required neutrino mass 
of 0.08 eV< <1.3 eV at the 68% C.L. (a < 0). In this case the predicted 
mass has a relatively strong dependence on the value of the universal radio 
background. Physically it is easy to understand the reason. The small radio 
background leads to a relatively large UHE7 fraction in the observations. They 
do not loose that much energy. Thus, smaller incoming UHEz/ energy and larger 
m-i, is needed to discribe the data. 

We performed a Monte Carlo analysis studying higher statistics. In the near 
future, the Pierre Auger Observatory will provide a ten times higher statistics, 
which reduces the error bars in the neutrino mass to about one third of their 
present values. 
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Fig. 6. Summary of the masses of the heaviest neutrino required in the Z-burst sce- 
nario, with their 1 a (solid) and 2 a (dotted) error bars, for the case of an extragalactic 
and a halo background of ordinary cosmic rays and for various assumptions about the 
diffuse extragalactic photon background in the radio band {a = 0,h — 0.71, Hm = 
0.3, = 0.7, Zmax = 2). From left: strong 7 attenuation, moderate and minimal uni- 
versal radio background. 



It should be stressed that, besides the neutrino mass, the UHEj/ flux at the 
resonance energy is one of the most robust predictions of the Z-burst scenario 
which can be verified or falsified in the near future. The required flux of ultrahigh 
energy cosmic neutrinos near the resonant energy should be detected in the near 
future by AMANDA, RICE, and the Pierre Auger Observatory, otherwise the Z- 
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burst scenario will be ruled out (cf. Fig. If such tremendous fluxes of ultrahigh 
energy neutrinos are indeed found, one has to deal with the challenge to explain 
their origin. It is fair to say, that at the moment no convincing astrophysical 
sources are known which meet the requirements for the Z-burst hypothesis, i.e. 
which have no or a negative cosmological evolution, accelerate protons at least 
up to 10^'^ eV, are opaque to primary nucleons and emit secondary photons only 
in the sub-MeV region. It is an interesting question whether such challenging 
conditions can be realized in BL Lac objects, a class of active galactic nuclei 
for which some evidence of zero or negative cosmological evolution has been 
found (see Ref. Q and references therein) and which were recently discussed 
as possible sources of the highest energy cosmic rays . 




10'8 10'9 1020 IQZ\ 1022 1023 

E [eV] 

Fig. 7. Neutrino fluxes, F — ^ 'l2i=ii-^"i required by the Z-burst hypothesis for 

the case of a halo and an extragalactic background of ordinary cosmic rays, respectively 
{a = 0,h = 0.71, Om = 0.3, i^A = 0.7, Zmax = 2). Shown are the necessary fluxes 
obtained from the fit results for the case of a strong UHE7 attenuation. The horizontal 
errors indicate the 1 a (solid) and 2 a (dotted) uncertainties of the mass determination 
and the vertical errors include also the uncertainty of the Hubble expansion rate. Also 
shown are upper limits from Fly's Eye on Fi,^ + Fv^ and the Goldstone lunar ultrahigh 
energy neutrino experiment GLUE on X^^^^. {F^^ + F^^), as well as projected sensi- 
tivities of AMANDA on F^^ + F^^ and Auger on F^^ + F^^ . The sensitiviy of RICE is 
comparable to the one of Auger. 
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